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Preface 


The  main  objective  of  this  effort  is  to  develop  a  code 
that  yields  a  reasonably  accurate  loading  diagram  for  a 
typical  high-work  gas  turbine.  A  practical  execution  time 
to  enable  a  designer  to  quickly  arrive  ac  an  acceptable 
blade  shape  is  a  secondary  objective. 

Two  people  deserve  special  mention  for  their  assistance 
and  encouragement.  I  am  particularly  indebted  to  Lt  Col 
Jack  Mattingly  of  the  Aero-Propulsion  Laboratory  for  sharing 
his  enthusiasm  and  knowledge.  It  was  a  joy  and  a  privilege 
to  be  his  student'.  The  success  of  this  project  is  largely 
due  to  Mr.  Robert  Gray,  the  project  sponsor,  also  of  the 
Aero-Propulsion  Laboratory.  Bob's  love  of  engineering, 
fascination  with  the  basic  physics,  and  mastery  of 
mathematics  made  this  effort  both  enjoyable  and  achievable. 

I  wish  to  thank  my  thesis  advisor.  Dr.  Ahmed  Halim, 
for  his  assistance  throughout  the  effort.  Thanks  also  to 
Dr.  Joseph  Shang  and  Lt  David  Amdahl  of  the  Flight  Dynamics 
Laboratory.  Dr.  Shang  was  willing  to  share  his  expertise 
and  provide  assistance  whenever  called  upon.  Lt  Amdahl 
graciously  provided  the  grid  generation  code  along  with 
instructions  for  its  use.  Finally,  the  greatest  debt  is 
owed  to  my  v^^£6l||||||||||||||||^for  her  patience  and  understanding 
through  the  many  days  and  nights  I  spent  in  front  of  a 
computer  terminal. 

Mark  A. Driver 
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Abstract 

A  numerical  algorithm  is  developed  with  the  capability 
of  capturing  shocks  in  the  internal  blade  passages  of  a 
modern  gas  turbine.  The  algorithm  uses  MacCormack's 
explicit  finite  difference  scheme  to  solve  the  two- 
dimensional  form  of  the  Euler  equations.  Inlet  and  exit 
boundary  conditions  are  developed  that  allow  disturbances  to 
propagate  out  of  the  computational  domain  without 
reflection.  Periodic  boundary  conditions  are  applied  such 
that  an  infinite  cascade  is  modeled. 

The  computed  steady  state  solution  is  compared  with 
experimental  data  for  a  high-work  low  aspect  ratio  turbine. 
The  ability  to  obtain  a  reasonably  accurate  blade  loading 
diagram  within  a  practical  execution  time  is  demonstrated. 
Two  oblique  shocks,  typical  of  those  formed  at  the  trailing 
edge  of  a  transonic  rotor  blade,  are  captured.  These  shocks 
are  smeared  over  several  grid  points,  as  expected  with  a 
shock  capturing  scheme,  but  their  influence  on  the  blade 
loading  diagram  Is  evident. 
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DEVELOPMENT  OF  A  SHOCK  CAPTURING  CODE  FOR  USE  AS  A  TOOL  IN 


DESIGNING  HIGH-WORK  LOW  ASPECT  RATIO  TURBINES 

Introduction 

Modern  aircraft  gas  turbine  engines  have  higher  core 
engine  pressures  and  temperatures  than  their  predecessors. 
These  engines  require  high  pressure  turbines  with  higher 
work  output  and  lower  corrected  mass  flow  than  earlier 
generation  turbines  (20:1).  This  type  turbine  incorporates 
blading  with  high  turning,  high  hub-to-tip  ratio,  and  low 
aspect  ratio.  The  high  hub-to-tip  ratio  and  low  aspect 
ratio  are  both  detrimental  to  turbine  efficiency.  The  high 
hub-to-tip  ratio  leads  to  boundary  layer  growth  on  the  inner 
and  outer  annulus  walls,  while  the  low  aspect  ratio  leads  to 
secondary  flow  losses  from  turning  the  boundary  layer 
through  a  large  angle  (4:267,271). 

The  flow  field  typically  consists  of  regions  of 
subsonic,  transonic,  and/or  supersonic  flow;  with  shock 
waves  often  forming  in  the  transonic  and  supersonic  regions. 
The  shock  waves  further  reduce  turbine  efficiency  due  to 
total  pressure  losses  and  shock-boundary  layer  interaction. 
Figure  1  is  a  schlieren  photograph  of  the  flow  through  a 
transonic  turbine  rotor  cascade  with  an  isentropic  exit  Mach 
number  of  1.15.  The  shocks  emanating  from  the  trailing 
edge  of  the  blades  and  the  reflection  from  the  suction 
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surface  are  clearly  visible,  as  is  the  presence  of  a 
separation  bubble  where  the  shock  intersects  the  suction 
surface.  This  figure  also  depicts  the  spatial  periodicity 
of  the  flow  field.  The  present  work  is  directed  toward 
locating  these  shocks  and  providing  the  designer  with  a 
reasonably  accurate  blade  loading  diagram. 

A  two-dimensional  solution  of  the  time  dependent  Euler 
equations  is  undertaken  since  they  provide  the  capability  to 
compute  mixed  subsonic,  transonic,  and/or  supersonic  flows 
and  also  capture  shocks.  A  modified  version  of  Shang's 
implementation  of  the  MacCormack  scheme  is  used  to  march  the 
solution  in  the  time  domain  (19). 

All  computations  are  carried  out  in  a  reference  system 
IB  attached  to  the  rotor.  Figure  2  depicts  a  typical  turbine 

stage  with  the  station  numbering  conventions  outlined  in 
reference  12.  The  subscript  R  denotes  the  relative 
reference  system,  defined  as  the  reference  system  moving 
with  the  rotor.  Velocities  in  this  system  are  given  the 
symbol  W,  while  velocities  in  the  nonrotating  reference 
system  are  given  the  symbol  V.  The  symbol  U  denotes  the 
rotational  speed.  A  computational  domain  is  established 
between  two  adjacent  rotor  blades.  The  domain  extends 
upstream  of  station  2^  and  downstream  of  station  3^  by 
approximately  one-half  of  the  rotor  axial  chord  length.  A 
typical  grid,  with  76  points  in  the  axial  direction  and  33 
points  in  the  tangential  direction,  is  shown  in  Figure  3. 
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The  code  developed  during  this  effort,  hereafter 
referred  to  as  BL02BLD,  is  used  to  investigate  three 
different  problems  with  known  steady  state  solutions.  A 
cascade  of  wedges  with  an  inlet  Mach  number  of  2.0  and 
completely  supersonic  flow  is  used  to  demonstrate  the 
ability  of  BLD2BLD  to  capture  well  defined  oblique  shocks. 
The  results  are  compared  to  the  exact  solution  presented  by 
Denton  (6:7).  For  the  second  case,  BLD2BLD  is  used  to 
compute  the  flow  field  in  the  rotor  passage  of  a  high-work 
low  aspect  ratio  turbine  tested  in  the  NASA  Lewis  Research 
Center's  Warm  Core  Turbine  Test  Facility  (20).  The  loading 
diagram  obtained  from  BLD2BLD  is  compared  to  the  design 
loading  diagram  obtained  from  NASA's  TSONIC  code  (20).  The 
C*  final  case  utilizes  the  NASA  turbine's  geometry  but  the  exit 

static  pressure  is  reduced  beyond  that  tested  by  NASA.  This 
case  is  thought  to  be  representative  of  the  conditions 
currently  being  investigated  by  turbine  designers. 
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Figure  3.  Typical  Grid  Used  In  the  Present  Analysis 
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II.  Analysis 


The  present  section  provides  a  synopsis  of  the 
governing  equations  and  boundary  conditions  applied  in  the 
BLD2BLD  code.  The  Euler  equations  are  written  in 
conservative  form  and  their  applicability  to  shock  capturing 
is  discussed.  A  coordinate  transformation  is  then  applied 
to  facilitate  the  numerical  solution  and  application  of 
boundary  conditions.  Next,  the  appropriate  boundary 
conditions  for  flow  through  a  blade  passage  are  developed. 
These  boundaries  are  the  radiative  boundaries  at  the  inlet 
and  exit,  the  periodic  fluid  boundaries,  and  the  blade 
surfaces . 

Governing  Equations 

The  Euler  equations  are  statements  of  the  conservation 
laws  for  mass,  momentum,  and  energy  assuming  an  inviscid 
nonconducting  gas.  The  favorable  pressure  gradient  and  high 
Reynold's  number  in  an  axial  flow  turbine  make  the  Euler 
equations  an  attractive  alternative  to  the  Navier-Stokes 
equations  as  long  as  heat  transfer  or  skin  friction  data  are 
not  desired.  The  customary  forms  of  the  conservation  laws 
for  mass,  momentum,  and  energy  in  the  absence  of  body  forces 
are 
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where  p  is  the  density,  W  is  the  velocity  vector,  p  is  the 
pressure,  e  is  the  internal  energy  per  unit  mass,  and  ^  is 
the  substantial  or  total  derivative  given  by 


D_ 

Dt 


at 


^  +  W*7 


(4) 


The  symbol  W  is  used  as  a  reminder  that  the  governing 
equations  are  written  in  a  frame  of  reference  attached  to 
the  rotor.  When  Eqs  (1),  (2),  and  (3)  are  expanded  using 
Eq  (4)  and  rearranged  such  that  p,  pu,  pv,  and  E^  are  the 
^4  dependent  variables  the  conservative  or  divergence  form  of 

the  governing  equations  is  obtained.  The  symbols  u  and  v 
denote  components  of  the  velocity  vector,  W,  in  the  x  and  y 
directions  defined  by  a  Cartesian  coordinate  system.  E^  is 
defined  to  be  the  total  energy  per  unit  volume  given  by 

E^  =  pj^e  +  u*/2  +  v*/2j  (5) 

Roache  states  that  the  conservative  form  of  the  governing 
equations  for  invlscid  compressible  flow  was  given  by 
Courant  and  Friedrichs  in  1946  (15:211).  Lax  showed  that 
the  conservative  form  of  the  governing  equations  satisfies 
the  weak  solution  of  the  Rankine-Hugonlot  relations  and  thus 
correctly  predicts  the  jump  conditions  across  the  shock 
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discontinuity  (15:211,  3:272).  In  fact,  use  of  the 
conservative  form  is  necessary  for  the  discontinuity  to 
physically  represent  a  shock  wave  when  numerical 
shock-capturing  schemes  are  applied  (3:272). 

References  7  and  15  present  derivations  of  the 
conservative  form  of  the  Euler  equations  and  such  a 
derivation  will  not  be  repeated  here.  The  results  will 
simply  be  stated  and  placed  in  vector  form  as  given  In 
reference  3. 

Eq  (1)  is  already  in  conservative  form  which  is 
apparent  when  Eq  (4)  is  used  to  expand  the  substantial 
derivative.  Some  manipulation  of  Eqs  (2)  and  (3)  is 
necessary  to  obtain  their  conservative  forms  with  the  end 
result  being,  for  two-dimensional  flow 


ItC^t) 


ItM  +  p)  + 

ItC^)  ^  p)  =  0  (8) 
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The  conservative  form  is  often  referred  to  as  the 
divergence  form  because  the  equations  identify  the 
divergence  of  physical  quantities  (3:50).  This  structure 
allows  the  governing  equations  to  be  written  in  vector  form: 
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where  U  contains  the  dependent  variables,  E  contains  the 
terms  differentiated  with  respect  to  x,  and  F  contains  the 
terms  differentiated  with  respect  to  y.  The  elements  of  V, 
E,  and  F  are: 


pu*+  p 


(Et+  p)u 


pv*+  p 
(E^+  p)v 


The  governing  equations  are  more  easily  solved  in  a 
rectangular  uniformly-spaced  computational  domain  than  in  a 
nonrectangular  physical  domain  with  nonuniform  spacing. 
Thus,  the  governing  equations  need  to  be  transformed  from  a 
Cartesian  coordinate  system  to  a  general  coordinate  system. 
The  general  spatial  transformation 


?  =  f(x,y) 
T?  =  T7(x,y) 


is  used  to  transform  Eg  (10)  from  the  physical  domain  (x,y) 
to  the  computational  domain  (^,n)>  The  partial  derivatives 
in  the  physical  domain  become 


^x  ^  +  ’^x  ^ 

3  9 


where  the  subscripts  x  any  y  denote  differentiation  with 
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respect  to  x  or  y.  The  terms  known 

X  y  X  y 

as  metrics.  The  Jacobian  of  the  transformation  given  by 
Eqs  (11)  and  (12)  is 


dTxry)  ■ 


X 

y| 

">x 

as 

follows : 

(15) 


J  =  1/J"*=  1/ 

where  the  subscripts  ?  and  t)  denote  differentiation  with 
respect  to  ?  or  t).  The  differentiations  are  carried  out 
using  finite  difference  representations  since  analytical 
expressions  are  not  available  (3:254).  The  metrics  are 
determined  from  the  relations 
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lyj  y,l  *  %  ^0 


(16) 
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(17) 

-J 

(18) 

y 

^x= 

-J  y^ 

(19) 

J  x^ 

(20) 

Applying  this  generalized  transformation  to  the  governing 
equations  in  vector  form,  Eq  (10),  results  in 

U.+  fE_  +  r)E+(rF,  +  77F=0  (2 

t  'x  ?  'x  77  ^y  ?  'y  r) 

A  numerical  solution  is  sought  for  the  governing  equations 
in  the  above  form. 
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Boundary  and  Initial  Conditions 


The  Euler  equations,  given  in  vector  form  in  Eq  (10), 
are  hyperbolic  for  all  flow  regimes  as  long  as  the 
time-dependent  terms  are  retained.  This  hyperbolic  behavior 
requires  that  both  boundary  and  initial  conditions  be 
specified.  The  importance  of  boundary  and  initial 
conditions  in  the  solution  of  partial  differential  equations 
(PDE's)  has  been  commented  on  by  Roache: 

A  first-order  ordinary  differential  equation  such 
as  df/dx  =  0  specifies  the  solution  of  a  problem  up  to 
an  additive  constant;  the  boundary  condition  determines 
the  value  of  the  constant .  A  first-order  partial 
differential  equation  such  as  df(x,y)/dx  =  0  specifies 
very  little  of  the  solution;  any  function  g(y) 
satisfies  the  PDE,  and  the  boundary  conditions  must 

specify  the  /unction,  A  PDE  such  as  =  C  Iw  is  the 

stream  function  and  C  is  the  vorticity}  really 
contains  very  little  information  on  y.  All  the 
fantastic  flow  patterns  of  common  gases  and  liquids  are 
solutions  of  the  sams  POE's,  the  Navier-Stolces 
equations.  The  flows  (solutions)  are  distinguished 
only  by  boundary  and  Initial  conditions,  and  by  the 
flow  parameters  such  as  Re  (Reynold's  number]  (15:139). 

Roach  suggests  that  the  most  difficult  boundary  condition  in 

compressible  flow  occurs  at  a  simple  wall  (15:261). 

However,  for  internal  flows  the  most  difficult  boundary 

conditions  occur  at  the  inlet  and  exit  of  the  passage  when 

the  flow  is  subsonic  at  one  or  both  of  these  locations. 

Also,  the  periodic  boundary  conditions  required  for  an 

infinite  cascade  model  are  a  formidable  challenge.  The 

boundaries  that  must  be  dealt  with  are  shown  in  Figure  4. 

Radiative  Boundary  Conditions.  The  current  effort  is 

directed  tovrards  achieving  a  steady  state  solution  in  the 
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Figure  4.  Physical  Boundaries 
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rotor  passage.  Unsteady  behavior  due  to  interaction  between 
the  stator  and  rotor  blade  rows  is  not  considered.  Since 
Eg  (10)  is  hyperbolic  it  exhibits  wave  behavior  (7:161). 

All  waves  radiating  outward  from  inside  the  computational 
domain  should  pass  the  inlet  or  exit  without  being 
reflected.  This  requirement  is  met  by  recasting  the 
governing  equations  into  characteristic  form  and  applying 
simple  wave  theory. 

Erdos  and  Alzner  (8:21-33)  and  Scott  and  Hankey 
(16:144-145)  report  use  of  the  method  of  characteristics  to 
propagate  disturbances  out  of  the  computational  domain. 
However/  Scott  and  Hankey  only  applied  the  method  at  the 
upstream  boundary.  They  dealt  with  a  compressor  cascade 
%f  having  a  supersonic  relative  velocity  entering  the  rotor 

with  the  axial  flow  remaining  subsonic.  The  relative 
velocity  exiting  the  rotor  was  subsonic.  They  inserted  a 
convergent -divergent  nozzle  downstream  of  the  exit  to 
achieve  supersonic  flow  so  that  supersonic  flow  conditions 
could  be  applied  at  the  outflow  boundary  (16:145).  The 
insertion  of  a  convergent-divergent  nozzle  downstream  of  the 
rotor,  while  a  novel  solution  for  the  compressor  problem, 
limits  the  exit  flow  from  the  rotor  to  subsonic  velocities. 
Modern  high-work  turbines  often  have  supersonic  exit 
velocities.  To  allow  generalization  of  the  current  work  to 
both  subsonic  and  supersonic  exit  velocities,  simple  wave 
theory  is  applied  herein. 
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Erdos  and  AIzner  state. 

The  characteristic  surfaces  formed  by  the 
hyperbolic  system  of  differential  equations  .  .  . 
consist  of  a  conoid  with  its  base  on  the  x,y  plane  and 
within  it  a  stream  path  which  intersects  the  conoid  at 
its  vertex.  If  the  vertex  is  placed  at  a  grid  point  at 
time  ti-At,  the  base  covers  the  domain  of  dependence  of 
the  point  at  time  t  (8:22). 

A  more  manageable  approximation  of  the  above  two-dimensional 
characteristic  theory  is  obtained  when  the  characteristic 
surfaces  are  replaced  by  characteristic  curves.  This  yields 
the  familiar  one-dimensional  characteristic  directions  shown 
in  Figure  5  (8:23).  The  inlet  and  exit  variables  that  are 
not  directly  specified  as  boundary  conditions  are  updated 
using  the  compatibility  relations,  or  Riemann  invariants, 
that  are  valid  along  the  appropriate  characteristics. 

The  characteristics  AO  and  CO  shown  in  Figure  5  are  wave 
paths  while  the  remaining  characteristic,  BO,  is  the 
particle  path.  The  wave  paths  are  given  by 


dx  . 

Jj-  =  u  .  a 

(22) 

^  .  u  -  a 

(23) 

where  a  is  the  local  speed  of  sound. 

The  remaining  equation 

is  that  of  the  particle  path: 

dx 

dt  ■  “ 

(24) 

Referencing  Erdos  and  AIzner  (8:22),  some  discussion  of 
the  characteristic  directions  is  in  order.  If  Station  2  is 
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Figure  5.  Characteristic  Directions 


the  Inlet  plane.  Station  3  lies  Inside  the  computational 
domain  and  Station  1  does  not  exist.  Only  one 
characteristic,  CO  In  Figure  6a,  lies  inside  the  domain  to 
relate  the  conditions  at  t”**  to  the  known  conditions  at  t”. 
Three  boundary  conditions,  described  later,  must  be 
specified  at  the  inlet.  One  of  these  boundary  conditions 
replaces  the  compatibility  relation  along  the  wave  path 
upstream  of  the  inlet.  The  other  two  conditions  replace  the 
characteristic  along  the  particle  path.  The  above 
statements  apply  when  the  Mach  number  at  the  inlet  is 
subsonic.  If  the  Mach  number  at  the  inlet  were  supersonic, 
all  characteristics  would  originate  upstream  of  the  inlet 


Figure  6.  i^ipllcable  Characteristics  at  Inlet  and  Exit 

m 

and  the  solution  at  this  boundary  would  be  known. 

Similarly,  if  Station  2  is  the  exit  plane.  Station  1 
lies  Inside  the  computational  domain  and  Station  3  does  not 
exist.  Two  characteristics,  AO  and  BO  in  Figure  6b,  lie 
inside  the  domain  to  relate  conditions  at  t'^^*  to  the  known 
conditions  at  t'*.  Only  one  boundary  condition  at  the  exit 
can  be  specified,  replacing  the  compatibility  relation  along 
the  wave  path  downstream. 

Particle  Path  Characteristics.  For  an  inviscld 
adiabatic  flow  in  the  absence  of  shocks,  an  alternate  form 
of  the  energy  equation  written  along  the  particle  path  is 
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where  s  is  the  entropy  (7:31).  Integration  of  Eq  (25)  along 
the  particle  path  given  by  Eq  (24)  yields  s  =  constant 
Thus  the  compatibility  relation  along  the  particle  path  can 
be  written  as 

p/p^~  constant  ( 26 ) 

Wave  Path  Characteristics.  The  derivation  of  the 
wave  path  characteristics  follows  the  outline  given  by  Erdos 
and  Alzner  (8:25-28).  For  clarity,  the  equation  of  state 
that  Erdos  and  Alzner  begin  with  will  be  derived.  This 
equation  of  state  relates  pressure  to  the  state  variables 
density  and  entropy.  The  starting  point  is  the  Gibbs 
equation: 

Tds  =  de  +  pdo  (27) 

where  T  is  the  temperature  and  v  is  the  specific  volume. 

The  static  enthalpy,  h,  is  by  definition 

h  =  e  +  pv  (28) 

thus 

de  +  pdv  =  dh  -  vdp  (29) 

Eq  (29)  is  substituted  into  Eq  (27)  giving 

Tds  =  dh  -  vdp  (30) 

Eq  (30)  is  rewritten  using  the  equation  of  state,  P  »  pRT  , 
and  the  differential  form  of  enthalpy  for  a  perfect  gas. 
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to  obtain 


dh  *  c  dT 


ds  =  C^dT/r  -  R  dp/p 


(31) 


where  is  the  specific  heat  at  constant  pressure  and  R  is 
the  gas  constant.  For  a  calorically  perfect  gas.  Eg  (31)  is 
integrated  between  states  1  and  2  to  give 


As/fe  =  (pyP*)j  <32) 

where  As  =  s^-  s^  and  y  is  the  ratio  of  specific  heats. 
Emanuel  then  perforins  the  substitution  (7:17) 


s^/R  =  constant 

leading  to 

p  =  kp^exp|s/cj 


(33) 


(34) 


where  )c  is  a  constant  and  the  state  2  subscript  has  been 
dropped.  Eg  (34)  is  the  state  eguation  used  to  initiate  the 
derivation  of  the  compatibility  relations  along  the  wave 
path  characteristics. 

Talcing  the  natural  logarithm  of  both  sides  of  Eg  (34) 
and  differentiating  the  result  using  the  substantial 
derivative  leads  to 

l/b  g|  =  y/o  gf  t  l/c,  H  (35) 


19 


where  C^ls  the  specific  heat  at  constant  volume.  Applying 
the  definition  of  the  substantial  derivative  and  using  the 
energy  equation,  Eq  (25),  gives 

i/t'[ ^ ^ ^ ]  ■  >'/p(  w* 


The  continuity  equation,  Eq  (6),  is  substituted  into  the 
above  relation,  and  the  result  multiplied  by  the  speed  of 
sound,  to  obtain 

The  next  step  is  to  add  and  subtract  the  streamwise  momentum 
equation.  Eg  (7),  to  the  above  relation.  Eq  (7)  is  first 
transformed  into  nonconservative  form  using  the  continuity 
equation: 


u 


£u 


+  V 


^  -  i/p  Si 


dy(. 


(38) 


Substituting  Eq  (38)  into  Eq  (37)  and  collecting  lilce  terms 


This  equation  becomes  an  exact  differential: 


a  dp  +  du  _  _  £v 
yp  fft  “  3F  “  9y 


(40) 
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on  the  equations  of  the  characteristic  lines  (8:26) 


dx  _  dy  _ 
u±a  “  V  “ 


dt 


(41) 


Eq  (41)  represents  the  characteristic  directions  shown  in 
Figure  5.  With  ds  =  0  ,Eq  (31)  yields  ^  =  pr  ^  • 

Substituting  this  relation  into  Bq  (40)  and  neglecting  the 
term  ^  outside  the  computational  domain  gives 


a  dT 
y-1  T 


du 


0 


(42) 


Eq  (42)  is  integrated  to  yield  the  well  known  Riemann 
invariants : 


^  ±  u  =  constant 

The  Riemann  invariants  are  often  referred  to  as  integrated 
forms  of  the  compatibility  relations  (7:163). 

Determination  of  Inlet  Boundary  Conditions.  The 
key  to  determining  the  proper  boundary  conditions  is  to 
treat  the  inlet  as  if  it  was  part  of  a  duct  extending 
infinitely  far  upstream.  All  waves  radiating  from  the 
computational  domain  should  pass  the  inlet,  without 
reflection,  and  continue  travelling  upstream  for  all  time. 

As  previously  mentioned,  if  Station  2  in  Figure  5  is  taken 
to  be  the  inlet  plane  then  the  only  straight  characteristic 
that  exists  is  the  line  CO.  This  straight  characteristic  is 
known  as  a  characteristic  of  the  first  kind  (9).  The 
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straight  characteristics  AO  and  BO  do  not  exist  since  the 
wave  is  travelling  upstream.  However,  a  curved 
characteristic  exists  which  extends  into  a  quiescent  region 
of  uniform  thermodynamic  properties.  This  characteristic  is 
known  as  a  characteristic  of  the  second  kind  (9).  These 
characteristics  intersect  at  point  O  in  Figure  6a,  a  grid 
point  at  the  inlet,  at  time  Application  of  the 

Riemann  invariants  along  both  types  of  characteristics 
results  in  a  system  of  two  simultaneous  equations  to  be 
solved  for  the  variables  u  and  a.  Specification  of  the 
velocity  direction  at  the  inlet,  along  with  the  assumption 
of  one-dimensional  flow  upstream  and  immediately  downstream 
of  the  inlet,  allows  use  of  a  coordinate  system  with  axes 
aligned  parallel  and  perpendicular  to  the  flow.  In  this 
coordinate  system  the  Riemann  invariants  given  by  Bq  (43) 
become 

±  W  =  constant  (44) 

where  W  is  the  magnitude  of  the  velocity  vector  W.  The 
benefits  of  the  one -dimensional  flow  assumption  immediately 
downstream  of  the  inlet  will  become  apparent  when  the 
numerical  implementation  of  the  characteristic  theory  is 
discussed  in  the  next  chapter.  The  two  simultaneous 
equations  are  now  in  the  variables  W  and  a. 

The  procedure  for  determining  the  inlet  boundary 
conditions  starts  with  specification  of  a,  W,  and  p  in  the 


quiescent  reqlon  upstream.  These  values  are  denoted  by  the 


symbols  a^  and  p^.  The  inlet  speed  of  sound  and  I 

velocity  at  time  t''**  are  related  to  the  known  speed  of 
sound  and  velocity  at  point  C  at  time  t'*  through  the  Riemann 
invariant  along  the  upstream  travelling  wave:  I 


2 

r-1 


a 


V  = 


2 

y-1 


V 


c 


(45) 


where  the  subscript  c  denotes  the  values  at  point  C. 
Along  the  characteristic  of  the  second  kind 


a  +  W  =  a  +  W  (46 

y-1  00  00  y-1 

Simultaneous  solution  of  Bqs  (45)  and  (46)  for  the  speed  of 
sound  at  the  inlet  results  in 

•  ■  H  V  ^  ( V  "c ) 


With  the  speed  of  sound  known,  the  velocity  is  obtained  from 
Bq  (46): 


FT  « 

The  pressure  is  obtained  from  the  isentropic  relation 

P  -  P»(  << 

With  the  speed  of  sound  and  pressure  known  the  state  point 
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Is  fixed,  uniquely  determining  the  density  and  internal 
energy. 


Determination  of  Exit  Boundary  Conditions.  The 
development  of  the  exit  boundary  conditions  closely  follows 
that  for  the  inlet  conditions.  The  exit  is  treated  as  an 
open-end  duct  that  exhausts  into  a  plenum.  This  requires 
that  the  exit  pressure  match  the  plenum  pressure.  Again 
referring  to  Figure  5,  the  characteristics  AO  and  BO 
originate  inside  the  domain  and  the  characteristic  CO  does 
not  exist  since  the  wave  is  travelling  downstream.  Since 
entropy  is  convected  along  particle  paths  (7:24),  the 
density  is  determined  from  the  isentropic  relation 

»b  (■  P/Pb  3*^’’  (50) 

where  p  has  been  specified  and  the  subscript  b  denotes  the 
values  at  point  B.  The  pressure  and  density  fix  the  state 
point  so  the  temperature,  hence  speed  of  sound  and  internal 
energy,  is  uniquely  determined.  The  exit  velocity  is 
obtained  by  evaluating  the  Riemann  invariant  along  the 
characteristic  AO; 


Lateral  Boundary  Conditions.  Only  one  blade  passage  of 
an  infinite  cascade  is  being  analyzed,  therefore  the 
conditions  along  the  lateral  boundaries  must  be  periodic 
except  on  the  blade  surfaces.  Along  any  axis  perpendicular 


to  the  axial  direction 


(pv),-  (pv)j^ 

where  1  denotes  the  lower  boundary  and  JL  denotes  the  upper 
boundary.  For  an  inviscid  flow,  the  only  condition  that  can 
be  applied  along  the  blade  surface  is  the  requirement  of 
surface  tangency. 

Initial  Conditions.  Two  classes  of  initial  conditions 
are  investigated.  The  first  class  is  termed  a  cascade 
tunnel  start  because  of  its  analogy  to  the  starting  of  an 
indraft  cascade  tunnel.  The  domain  is  initialized  at  the 
stagnation  conditions  associated  with  the  quiescent  region 
upstream  of  the  inlet.  At  time  t^  the  exit  pressure  is 
instantaneously  applied  just  as  if  a  valve  had  been  opened 
to  a  low  pressure  plenum  downstream  of  the  exit.  The  second 
class  of  initial  conditions  consists  of  using  a  restart  file 
available  from  a  previous  solution.  The  restart  file  is 
used  to  initialize  the  domain  before  applying  new  quiescent 
region  conditions  or  a  new  exit  pressure.  As  an  example,  if 
solutions  are  desired  for  several  different  values  of  the 
exit  pressure,  the  cascade  tunnel  start  could  be  used  for 
the  highest  exit  pressure  and  then  the  restart  file  used  to 


(52) 

(53) 

(54) 

(55) 
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Initialize  the  domain  for  the  next  pressure  investigated,  or 
vice  versa. 

Of  the  two  classes  of  initial  conditions,  the  cascade 
tunnel  start  is  perhaps  the  most  noteworthy.  It  has  been 
pointed  out  by  Gray  (9)  that  the  cascade  tunnel  start  is 
possible  due  to  application  of  characteristic  theory  at  the 
radiative  boundaries.  If  the  dependent  variables  were 
explicitly  stated  at  the  inlet  and  exit  a  better  initial 
approximation  of  the  steady  state  solution  would  be  required 
such  that  large  disturbances  could  not  pass  these  radiative 
boundaries . 
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III.  Solution  Technique 


The  time-dependent  Euler  equations,  given  in  vector 
form  as  Eq  (21),  are  hyperbolic  with  respect  to  time.  This 
hyperbolic  nature  alloi/s  the  solution  to  be  marched  in  time 
even  though  the  flow  may  consist  of  mixed  subsonic, 
transonic,  and/or  supersonic  regions  (3:259).  The  numerical 
technique  used  to  obtain  a  solution  of  the  Euler  equations 
for  steady  flow  in  a  turbine  cascade  is  discussed  in  this 
chapter.  Details  of  the  finite  differencing  scheme  and 
numerical  implementation  of  the  boundary  conditions  are 
covered.  Stability  and  convergence,  allowable  time  step 
size,  and  numerical  damping  are  also  mentioned. 

Finite  Difference  Scheme 

MacCormack's  explicit  finite  difference  scheme  is  used 
to  obtain  a  time  marching  solution  of  the  Euler  equations 
(21).  The  BLD2BLD  code  is  a  modified  version  of  Shang's 
three-dimensional  code  (17,18).  The  modifications  were 
necessary  to  adapt  the  code  for  solution  of  two-dimensional 
internal  flows.  These  modifications  have  no  effects  on 
Shang's  implementation  of  the  differencing  scheme  except  at 
the  lateral  boundaries.  MacCormack's  scheme  is  a  two-step 
Lax-Wendroff  method.  This  scheme  uses  alternating  forward 
and  backward  differences  for  the  predictor  and  corrector 
sweeps  respectively.  predictor  and  corrector  sweeps 

applied  to  the  Euler  equations  in  a  Cartesian  coordinate 


system,  Bq  (12),  are  defined  as  (11:151) 


I 


I 


I 


I 


U""*  =  u^'  -  -  E"  1  -  -  F''  1 

V,J  V.J  Ax  '••jJ  Ay^  v.j+l.  v.jJ 


(56) 


= 


1 


i 


(57) 


where  and  equal  E(u”^  and  F(u!^p.  The 
subscripts  i  and  j  denote  a  mesh  of  points  (x^,yp  with 
uniform  spacing  Ax  and  Ay.  The  superscript  n  refers  to  the 
time  t  s  nAt  where  At  is  the  time  increment  for  one 
predictor-corrector  cycle.  Finally,  the  overbar  denotes 
predicted  quantities.  MacCormack  points  out  that  using  a 
forward  difference  for  the  predictor  and  a  backward 
difference  for  the  corrector  represents  only  one  of  four 
possible  differencing  methods  (11:151).  In  this  effort  only 
the  variation  given  by  Eqs  (56)  and  (57)  is  implemented. 

The  code  achieves  maximum  efficiency  by  utilizing  the 
Courant-Friedr ich-Lewy  (CFL)  condition  derived  by  Shang 
from  a  stability  analysis  (18:4): 


(58) 
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where  At 


is  the  allowable  time  step  size  and  and  u 
crt  K  '0 

are  the  contravar lant  velocities  defined  as 

+  ?yV  (59) 

=  TJjjU  +  T)yV  (60) 

The  CFL  condition  requires  that  the  physical  domain  of 
dependence  lie  within  the  numerical  domain  of  influence. 

The  numerical  domain  can  be  larger,  but  not  smaller,  than 
the  physical  domain  (3:76).  Anderson  reports  that 
hyperbolic  systems  tend  to  behave  best  when  the  numerical 
and  physical  domains  of  influence  are  nearly  equal  (3:77). 
The  CFL  numbers  used  in  the  present  wor)c  ranged  from  0.2  to 
0.8,  with  0.8  used  most  often.  Satisfaction  of  the  CFL 
condition  is  necessary  for  stability  of  the  MacCormacic 
scheme.  Lax's  equivalence  theorem  equates  stability  and 
convergence  of  a  consistent  finite  difference  scheme,  for  a 
linear  system  of  equations  (3:49).  This  theorem  has  never 
been  proven  for  nonlinear  equations,  such  as  the  Euler 
equations,  but  this  work  proceeds  under  the  assumption  that 
the  theorem  applies. 

Explicit  numerical  damping  is  applied  to  suppress 
numerical  oscillations.  The  damping  is  applied  in  each 
sweep  direction  using  Shang's  modification  of  MacCormack's 
fourth-order  pressure  damping  terms  (17:1349): 

P^t  ♦  [?*  *  «y  r**]  P 
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where  ft  is  the  damping  coefficient.  In  the  present  work  a 
damping  coefficient  of  2  is  consistently  used. 

Implementation  of  Boundary  Conditions 

Attention  is  now  focused  on  the  numerical 
implementation  of  the  boundary  conditions  developed  in  the 
previous  chapter.  The  only  boundaries  that  require 
additional  computations  beyond  the  predictor -corrector 
differencing  are  the  inlet  and  exit.  The  conditions  at  the 
lateral  boundaries,  both  the  periodic  boundaries  and  the 
blade  surfaces,  are  obtained  with  modifications  of  the 
I  ^  predictor-corrector  differencing  scheme. 

Radiative  Boundaries.  The  simple  wave  theory  used  to 
determine  the  appropriate  values  of  the  dependent  variables 
I  at  the  inlet  and  exit  was  developed  in  Chapter  II.  The 

numerical  Implementation  of  this  theory  is  the  subject  of 
this  section.  Referring  again  to  Figure  6a,  if  station  2  is 
the  inlet  plane  the  conditions  at  any  node  point  at  time 
t”"*^  is  related  to  conditions  at  point  C  at  time  t”  through 
the  Riemann  invariants,  Eq  (45).  The  conditions  at  any 
inlet  node  point  are  similarly  related  to  conditions  in  the 
quiescent  region  upstream  through  Eq  (46).  The  speed  of 
sound,  velocity,  and  static  pressure  in  the  quiescent  region 
are  specified  as  inputs  to  the  code.  The  conditions  at 


*  (  "x  *  ’’H  P 
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point  C  thus  need  to  be  determined  in  order  to  evaluate  the 
conditions  at  the  inlet.  This  is  the  point  where  the 
one-dimensional  flow  assumption  immediately  downstream  of 
the  inlet  becomes  advantageous.  A  coordinate  system  is 
aligned  with  the  flow  direction  at  the  inlet  such  that  one 
axis,  denoted  X',  is  aligned  with  the  flow.  By  construction 
this  axis  is  also  aligned  along  a  line  of  r),  hence  the  j 
index,  equal  to  a  constant  as  shown  in  Figure  3.  This 
construction  allows  the  equation  of  the  characteristic, 

Eq  (23),  to  be  written  as 

dx'/dt  =  W  -  a  (63) 

G'  C  G 

Point  x^is  located  through  an  iterative  solution  of  the 
equation  obtained  by  integrating  Bq  (63): 

x'  =  xr  ,  -  r  W  -  a  )At  (64) 

with  linear  interpolation  used  to  determine  the  flow 
properties  at  point  C.  The  code  uses  the  Newton-Raphson 
method  (10:764)  to  perform  the  iterative  solution  of  Eq  (64). 
In  addition,  bisection  is  used  as  a  fall-safe  whenever  the 
Ne%fton-Raphson  method  would  take  the  solution  out  of  bounds 
(14:258).  Once  point  C  is  located  and  the  flow  properties 
obtained,  Eqs  (47)-(49)  are  used  to  determine  the  flow 
properties  at  the  inlet. 

The  exit  boundary  conditions  are  implemented  in  a 
similar  fashion  to  those  at  the  inlet.  Referring  to 
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Figure  6b,  If  station  2  is  taken  as  the  exit  plane  the 
characteristics  AO  and  BO  are  available  to  relate  conditions 


c« 


I 


I 


inside  the  domain  at  time  t'^  to  exit  conditions  at  time 
t"^*.  After  applying  a  one -dimensional  flow  assumption 
immediately  upstream  of  the  exit  the  characteristic 
equations  become 


c'  =  x'  -  C  W  +  a  1 

A  AXkl  ^  A  A 


At 


and 


x'  =  x'  .  -  W  At 

B  •xtl  B 


(65) 


(66) 


Points  A  and  B  are  located  through  an  iterative  solution  of 
Eqs  (65)  and  (66)  with  linear  interpolation  used  to 
determine  the  flow  properties.  With  the  exit  pressure,  P  , 
given  as  input  to  BLD2BLD  the  exit  boundary  conditions  are 
determined  using  Eqs  (50)  and  (51). 

Lateral  Boundaries.  The  flow  properties  along  the 
lateral  boundaries  are  not  explicitly  specified  but  rather 
are  determined  using  slightly  modified  forms  of  the 
MacCormack  scheme.  Along  the  periodic  boundaries  the 
conditions  given  by  Eqs  (52)-(55)  are  met  by  actually 
performing  the  predictor  and  corrector  sweeps  at  the 
boundary.  If  a  forward  predictor  is  applied  at  the  upper 
boundary,  denoted  by  JL,  the  flow  properties  at  j  =  JL  +  1 
are  needed.  Since  this  line  of  points  does  not  exist  when 
only  one  blade  passage  of  an  infinite  cascade  is  being 
modeled,  the  properties  at  j  »  jl  -i-  1  are  replaced  by  the 


32 


conditions  at  j  »  2  .  An  analogous  procedure  is  followed 
when  a  backward  predictor  is  applied  at  the  lower  boundary. 
If  the  lower  boundary  is  denoted  by  j  =  1  ,  the  conditions 
at  j  -  0  are  replaced  by  the  condition  at  j  =  JL  -  1  . 
This  differencing  procedure  appears  to  satisfy  the 
periodicity  requirement  up  to  approximately  one  point  away 
from  the  leading  and  trailing  edges. 

The  appropriate  boundary  condition  along  the  blade 
surfaces  is  the  surface  tangency  requirement,  W  •  n  =  0  , 
where  n  is  the  surface  normal  (15:261).  The  MacCormack 
scheme  is  applied  directly  at  the  blade  surface,  using 
forward  differences  for  both  the  predictor  and  corrector  on 
the  suction  surface  and  backward  differences  for  both  sweeps 
along  the  pressure  surface.  Differencing  in  the  same 
direction  for  both  the  predictor  and  corrector  is  usually 
unstable,  but  Anderson  points  out  that  the  surface  tangency 
condition  affects  the  solution  such  that  stability  is 
maintained  (3:277). 

Solutions  to  a  variety  of  test  cases  were  attempted 
prior  to  attempting  the  solution  for  an  infinite  cascade. 
These  test  cases  ranged  from  fully  subsonic  flow  through  a 
curved  duct  to  fully  supersonic  flow  through  a  duct  with 
surface  discontinuities  to  generate  shock  and  expansion 
waves.  The  above  boundary  conditions,  less  the 
periodicity  constraint,  were  applied  to  the  test  cases  with 


excellent  results. 


Numerical  difficulties  are  encountered  at  the  periodic 
downstream  boundaries  if  the  cascade  tunnel  start  is  used. 

A  similar  difficulty  was  observed  by  Scott  and  Hankey 
(16:146).  They  alleviated  the  problem  by  treating  the 
do%mstream  lateral  boundaries  as  solid  surfaces  until  the 
solution  evolved  to  a  point  that  the  flow  became  aligned 
with  the  channel  (16:146).  The  same  remedy  is  applied  to 
the  current  effort,  eliminating  the  difficulties. 


j 
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IV.  Results  and  Conclusions 

Results 

This  section  describes  the  results  obtained  when  the 
BLD2BLD  code  Is  applied  to  three  different  test  cases.  The 
first  solution  is  for  a  cascade  of  wedges  with  an  inlet  Mach 
number  of  2.0  and  supersonic  flow  throughout  the  passage.  A 
solution  is  then  obtained  for  the  NASA  high-work  low  aspect 
ratio  turbine  rotor  described  in  reference  20.  The  NASA 
turbine  geometry  is  used  for  the  third  solution  but  the  exit 
pressure  is  dropped  below  that  of  the  NASA  tests.  This 
results  in  a  supersonic  exit  Mach  number  typical  of  that 
currently  under  consideration  by  turbine  designers. 

Case  One.  Cascade  of  Wedges.  The  cascade  of  wedges 
presented  by  Denton  is  used  to  demonstrate  the  ability  of 
BLD2BLD  to  capture  well  defined  oblique  shocks  (6:7).  This 
cascade  is  shown  in  Figure  7.  The  cascade  has  an  inlet  Mach 
number  of  2.0  and  is  designed  such  that  the  leading  edge 
shock  is  exactly  canceled  upon  reflection  to  the  upstream 
corner  resulting  in  uniform  flow  between  the  two  parallel 
surfaces.  The  grid  used  consists  of  125  points  in  the  « 

axial  direction  and  26  points  in  the  tangential  direction. 
The  computed  static  pressure  contours  are  shown  in  Figure  8. 
The  shock  that  forms  at  the  leading  edge  is  smeared  over 
several  points,  as  is  typical  for  a  shock  capturing 
scheme,  and  is  smeared  even  more  upon  reflection.  As  a 


Figure  7.  Cascade  of  Wedges 


Static  Pressure  Contours  for  Cascade  of  Wedges 


result,  the  reflected  shock  does  not  exactly  meet  the 
upstream  corner  and  a  weak  shock  is  reflected  back  across 
the  passage.  Additionally,  a  weak  expansion  emanates  into 
the  region  between  the  parallel  surfaces.  Two  oblique  shock 
waves  occur  at  the  trailing  edge,  similar  to  the  oblique 
shock  waves  occurring  at  the  trailing  edge  of  the  turbine 
blades  shown  in  Figure  1.  These  trailing  edge  shocks  are  a 
result  of  the  periodicity  conditions  applied  immediately 
downstream  of  the  trailing  edge.  Figure  9  contrasts  the 
exact  and  computed  solutions  in  terms  of  Mach  number  versus 
the  non-dimensional  chord  length.  Numerical  oscillations, 
or  dispersive  errors,  occur  nt^ar  the  points  where  the  shock 
and  expansion  waves  are  generated  or  reflected.  These 
dispersive  errors  are  typical  of  second-order  methods  such 
as  HacCormack's  (3:92).  The  oscillations  near  the  leading 
edge  may  also  be  aggravated  by  the  periodicity  conditions 
immediately  upstream.  Figure  9  also  shows  that  the  weak 
shock  generated  due  to  smearing  intersects  the  lower  surface 
near  60  percent  of  the  chord  length.  An  exact  solution  for 
the  expansion  along  the  lower  surface  was  not  presented  by 
Denton.  Denton  attributes  the  exact  solution  to  Brown 
Boveri  &  Co.  of  Baden,  Switzerland  (6:9). 

Case  Two.  NASA  High-Work  Low  Aspect  Ratio  Turbine. 

The  BLD2BLD  code  is  now  used  to  obtain  a  steady  state 
solution  for  the  flow  in  a  blade  passage  of  an  experimental 
turbine.  The  experimental  turbine  is  a  0.767  scale  model  of 
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the  first  stage  of  a  two  stage  high  pressure  turbine 
designed  for  use  in  a  high  bypass  ratio  engine.  This  model 
was  tested  in  the  NASA  Lewis  Research  Center's  Warm  Core 
Turbine  Test  Facility  and  the  results  are  well  documented  in 
reference  20.  Even  though  strong  secondary  flo%rs 
occurred  in  the  tests  (20:5),  this  turbine  was  chosen  for 
comparison  because  of  the  difficulty  in  obtaining 
two-dimensional  cascade  data. 

Figure  10  shows  the  mean-line  velocity  diagram  obtained 
from  the  NASA  experiment.  Using  the  mean-line  blade 
coordinates  given  in  Table  1,  and  the  relative  gas  angles 
from  Figure  10,  the  grid  shovm  in  Figure  3  is  constructed. 
The  grid  consists  of  76  points  in  the  axial  direction 


Figure  10.  Experimental  Mean-Line  Velocity  Diagram  (20:22) 
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and  33  points  in  the  tangential  direction.  The  grid 
extends  approximately  one-half  axial  chord  length  upstream 
and  downstream  of  the  blades,  and  the  blade  pitch  is 
3.059  cm.  The  lines  of  constant  n  are  aligned  with  the 
relative  gas  angles  at  each  end  of  the  grid.  The  blade 
surface  coordinates  are  obtained  using  a  cubic  spline 
interpolation  (14:86-89).  In  addition,  the  rounded  trailing 
edge  was  replaced  by  a  sharp  trailing  edge  for  computational 
purposes.  The  grid  was  generated  using  the  algebraic 
portion  of  Amdahl's  ORTHGNL  grid  generation  code  (1). 

The  inlet  and  exit  conditions  necessary  for  input  to 
BLD2BLD  are  derived  from  the  following  NASA  test 
data  (20:6-9,  22): 


r  *  1.4 

P.  =  31.03  X  10*  N/ta* 
^1 

P.  /P,  =  1.652 
'2R 

C''Acr)2  “0-888 
C“/''cr)2  ■ 


T.  =  422.2  K 

ti 

P  /P2  »  1.704 


P.  /P.  =  2.360 

Cl  Cj 


(vAcr)3  =  0*3®^ 
(wAer)3  = 


where  t  identifies  total  properties,  R  denotes  a  frame  of 
reference  moving  with  the  rotor,  cr  is  a  condition 
corresponding  to  a  Mach  number  of  unity,  and  the  station 
numbering  and  velocity  symbols  correspond  to  those  shown  in 
Figure  2.  Note  that  the  test  data  specifies  velocity  ratios 
rather  than  Mach  numbers.  This  is  a  convention  typically 
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followed  by  turbine  designers.  Although  assumed  constant 
for  the  present  work,  the  specific  heats  can  vary 
significantly  through  some  turbines.  This  results  in 
significant  y,  hence  Mach  number,  variations.  Thus  velocity 
ratios  are  more  meaningful  to  the  designer.  The  Mach  number 
is  related  to  the  velocity  ratio  through  the  equation  (2:5) 


I 


I 


As  stated  in  Chapter  III,  the  speed  of  sound,  velocity,  and 
static  pressure  in  the  quiescent  region  upstream  of  the 
inlet  are  required  as  input  to  BLD2BLD.  The  speed  of  sound, 
velocity,  and  static  pressure  at  station  2R,  obtained  from 
the  test  data,  are  input  as  the  conditions  existing  in  the 
quiescent  region.  Assuming  the  flow  through  the  stator  is 
adiabatic,  the  static  temperature  at  station  2R  is  obtained 
from  (2:5) 


(70) 


resulting  in  T_  =  T-_  =  366.71  K  and  a_  =  383.91  m/s  . 

Bq  (69),  with  replaced  by  is  used  to  arrive  at 

**2R  ”  0*352  thus  giving  =  M2jj  =  135.14  m/s  .  The 
static  pressure  is  obtained  directly  from  the  given  values 
of  and  P^.  /p2  as  P^  =  18.21  X  10*  N/m*  .  The  static 
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pressure  at  station  3R  Is  input  to  BLD2BLD  as  the  exit 
pressure,  P^.  The  static  pressure  is  obtained  from  the 
relation  (2:5) 


P.  is  obtained  directly  from  the  test  data  giving 

Pg  =  Pjjj  =  12.05  X  10*  N/in*  .  All  required  input  to  BLD2BLD 
is  now  available. 

The  BLD2BLD  solution  shows  that  the  rotor  is  choked  at 
the  throat,  as  shown  by  the  Mach  contours  of  Figure  11  and 
the  static  pressure  contours  of  Figure  12.  The  NASA  test 
data  also  suggest  that  the  rotor  was  choked,  or  very  nearly 
so,  at  the  conditions  described  above  (20:3).  The  average 
value  of  (»Act)  2  is  slightly  lower  than  the  test  data, 

0.364  compared  to  0.381,  probably  due  to  a  slight  narrowing 
of  the  throat  dimension  when  the  blade  surfaces  were 
numerically  generated.  The  Mach  contours  of  Figure  11  and 
the  static  pressure  contours  of  Figure  12  suggest  that  two 
shock  waves  are  emanating  from  the  blade  trailing  edge.  The 
presence  of  these  shocks  is  verified  upon  examination  of 
Figures  13  and  14.  These  figures  show  that  the  shock  wave 
emanating  from  the  pressure  surface  strikes  the  suction 
surface  near  65  to  70  percent  of  the  axial  chord.  The  exact 
location  cannot  be  determined  because  of  shock  smearing.  A 
reflection  from  the  suction  surface,  similar  to  Figure  1, 
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TANGENTIAL  DISTANCE  (CM) 


TANGENTIAL  DISTANCE 


Figure  12.  Static  Press 


P/P2 


Figure  13.  Blade  Loading  Diagram  for  Case  Two 


x/c 


Figure  14.  Surface  Pressures  (P/P2)  ^or  Case  Two 


can  not  be  detected  because  of  this  smearing.  The  existence 
of  these  shocks  can  not  be  verified  from  the  experimental 
data,  since  NASA's  experimental  apparatus  was  not  equipped 
to  sense  the  presence  of  shocks.  Figures  11  and  12  also 
show  the  presence  of  a  second  shock  emanating  from  the 
trailing  edge  of  the  suction  surface.  These  shocks  form 
because  of  the  periodicity  requirement  downstream  of  the 
blade  trailing  edge.  If  the  lateral  boundaries  upstream  and 
downstream  of  the  blade  are  treated  as  solid  walls  these 
trailing  edge  shocks  are  not  generated.  Referring  to  Figure 
15,  a  normal  shock  then  forms  downstream  of  the  throat,  as 
verified  by  a  solution  using  BLD2BLD,  such  that  the 
specified  exit  pressure  is  met.  The  pressure  is 
discontinuous  across  the  normal  shock,  with  the  upstream  and 
downstream  pressures  denoted  by  p^^  and  p^^^  respectively. 

If  periodicity  is  then  applied,  a  pressure  discontinuity 
exists  across  the  stagnation  streamline.  The  normal  shock 
is  thus  replaced  by  two  oblique  shocks  that  align  themselves 
such  that  no  pressure  discontinuity  exists  across  the 
streamlines . 

Referring  again  to  Figure  13,  the  loading  diagram 
obtained  from  BLD2BLD  is  compared  with  the  design  loading 
diagram  obtained  from  NASA's  TSONIC  code.  The  two  codes 
compare  favorably  except  near  the  leading  edge  of  the 
pressure  surface  and  in  the  vicinity  of  the  shock  waves. 
Experimental  measurements  on  transonic  turbine  blades. 
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presented  by  Denton  (6:7),  suggest  that  the  BLD2BLD  solution 
is  more  representative  o£  the  blade  loading  near  the 
pressure  surface  leading  edge.  Comparisons  with  TSOMIC  can 
not  be  made  in  the  vicinity  of  the  shock  waves  for  two 
reasons.  First,  TSONIC  may  have  been  applied  to  a  geometry 
with  a  slightly  larger  throat  dimension,  thus  alleviating 
choking.  Second,  there  is  some  question  in  the  design 
community  about  TSOMIC's  ability  to  deal  with  shocks  (9).  A 
qualitative  comparison  of  BLD2BLD's  shock  capturing  ability 
can  be  made  against  the  widely  used  Denton  code  (6).  The 
Denton  code  employs  a  pitchwise  smoothing  of  the  flow 
properties  in  addition  to  pressure  damping.  This  causes  the 
shocks  to  be  smeared  to  the  extent  that  the  effect  of  a 
C*  shock  emanating  from  the  pressure  surface  is  almost  entirely 

damped  out  before  reaching  the  suction  surface.  Denton  has 
demonstrated  this  effect  even  at  exit  Mach  numbers  as  high 
as  1.42  (6:7).  Even  though  BLD2BLD  has  a  smearing  effect, 
as  do  all  shock  capturing  schemes,  the  effect  of  a  shock 
impinging  on  the  suction  surface  is  evident  from  Figures  13 
and  14.  Examination  of  Figure  16  reveals  a  reduction  in  the 
magnitude  of  the  velocity  vectors  near  the  suction  surface 
due  to  a  shock  extending  across  the  passage.  The  velocity 
magnitude  is  further  reduced  as  the  suction  surface  trailing 
edge  shock  is  encountered.  The  average  value  of 
from  the  BLD2BLD  solution  is  slightly  larger  than  the  value 
obtained  from  the  NASA  test  data,  0.879  using  BLD2BLD 
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Figure  16.  Velocity  Vectors  for  Case  Two 
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compared  to  an  experimental  value  of  0.841.  The  calculated 

exit  relative  total  pressure,  P.  ,  is  also  slightly  larger 

^3R 

4  2 

than  the  experimental  value,  19.53  X  10  M/n  calculated 
versus  18.64  X  10^  experimental.  This  is  not  surprising 

in  light  of  the  strong  secondary  flows  observed  in  the 
experiment . 

Case  Three.  NASA  Turbine  with  Lowered  Exit  Pressure. 
This  final  case  uses  the  case  two  geometry  and  inlet 
conditions,  but  the  exit  pressure  is  lowered  to 
8.05  X  10^  N/m^.  The  exit  pressure  is  chosen,  rather 
arbitrarily,  such  that  the  exit  Mach  number  is  supersonic. 
This  case  uses  the  same  grid  generated  for  case  two.  This 
case  is  thought  to  be  representative  of  the  conditions 
currently  under  consideration  by  turbine  designers. 

The  BLD2BLD  solution  yields  an  average  value  of 
("Act)  2  of  1.127  corresponding  to  an  exit  Mach  number  of 
1.159.  Figures  17  and  18  show  the  computed  Mach  and 
pressure  contours  respectively.  Shock  waves  are  again 
present  at  the  trailing  edges  of  the  pressure  and  suction 
surfaces.  Figures  19  and  20  show  that  the  shock  emanating 
from  the  pressure  surface  is  somewhat  stronger  than  for  case 
two.  The  shock  also  impinges  on  the  suction  surface  farther 
downstream  than  for  the  previous  case.  The  impingement  now 
occurs  at  approximately  80  percent  of  the  axial  chord.  The 
likelihood  of  flow  separation  increases  as  the  Impingement 
point  moves  toward  the  trailing  edge.  The  flow  tends  to 
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Figure  17.  Mach  Contours  for  Case  Three 
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Figure  20.  Surface  Pressures  (p/p^)  for  Case  Three 


reattach  when  the  shock  impinges  farther  upstream,  as  shovm 
in  Figure  1.  For  a  typical  turbine  rotor,  the  flow 
separates  from  the  suction  surface  just  prior  to  reaching 
the  trailing  edge.  Therefore  as  the  shock  impingement  point 
moves  closer  to  the  trailing  edge  the  possibility  of 
separation  increases.  Figure  21  shows  the  velocity  vectors 
for  case  three. 

Conclusions 

A  two-dimensional,  time  marching  Euler  code  has  been 
developed  to  obtain  steady  state  solutions  for  flows  through 
transonic  turbine  cascades.  The  simple  wave  theory  applied 
at  the  radiative  boundaries  allows  any  streamwise  travelling 
disturbance  to  propagate  out  of  the  domain  without 
reflection.  The  method  of  applying  the  periodic  boundary 
conditions  appears  to  correctly  model  the  conditions 
existing  in  an  infinite  cascade.  The  code  predicts  the 
formation  of  two  oblique  shocks  at  the  blade  trailing  edges, 
typical  of  those  experimentally  observed  in  transonic 
turbine  blades.  The  solutions  obtained  using  the  BLD2BLD 
code  compare  favorably  with  the  experimental  data  within 
certain  limits.  These  limits  become  evident  by  comparing 
computed  two-dimensional  data,  based  on  an  inviscid 
analysis,  with  three-dimensional  experimental  data  strongly 
Influenced  by  viscous  effects. 

The  data  processing  rate  is  2.6  X  10  ’  seconds  per  grid 
point  per  time  step  for  the  CRAY-XMP  computer.  The  solution 
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is  monitored  until  consecutive  calculations  show  less  than  a 


0.02  %  change  in  the  dependent  variables.  The  solution  is 
then  considered  to  be  the  asymptotic  steady  state  solution. 

A  typical  calculation  begun  with  the  cascade  tunnel  start 
Initial  conditions  regulres  approximately  8000  time  steps  to 
achieve  steady  state  convergence.  This  translates  to  8.69 
minutes  of  CRAY-XMP  CPU  time  for  the  76  X  33  grid  used  for 
the  NASA  turbine  case. 
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V.  Recommendations 


The  areas  suggested  for  further  study  can  be  divided 
into  two  categories,  those  related  to  improving  the  current 
ability  of  the  code  to  obtain  steady  state  blade-to-blade 
solutions  and  those  related  to  extending  the  capabilities  to 
handle  unsteady  flows  and  viscous  phenomena. 

Three  specific  recommendations  apply  to  the  first 
category.  They  are 

1.  Investigate  the  application  of  two-dimensional 
characteristic  theory  at  the  periodic  boundaries. 
This  should  allow  the  removal  of  disturbances 
propagating  in  the  tangential  direction  with 
minimal  application  of  explicit  numerical  damping. 

A  reduction  in  the  damping  may  allow  better  shock 
wave  definition  and  also  allow  capturing  of  the 
reflected  shock  from  the  suction  surface. 

2.  Remove  the  one-dimensional  flow  assumption  at  the 
inlet  and  exit  boundaries.  Characteristic  theory 
would  still  be  applied,  but  in  a  more  general 
fashion. 

3.  The  effect  of  grid  refinement  at  the  blade  leading 
and  trailing  edges  was  not  thoroughly  investigated 
during  this  effort.  Grid  refinement,  especially  at 
the  leading  edge,  should  allow  better  resolution  of 
the  flow  field  near  the  stagnation  points. 


The  recommendations  relating  to  the  second  category 
will  require  much  more  effort  than  those  in  the  first.  Two 


specific  areas  are  suggested  for  implementation: 

1.  Viscous  effects  need  to  be  accounted  for.  Low 
aspect  ratio  turbines  are  strongly  influenced  by 
viscous  phenomena  and  this  needs  to  be  accounted 
for  in  order  to  obtain  accurate  efficiencies. 

2.  Extend  the  code  to  solve  the  unsteady  vane-blade 
interaction  problem.  This  is  the  most  extensive 
effort  suggested  for  implementation.  Vane-blade 
interaction  is  currently  receiving  much  attention 
in  the  design  community  and  a  code  that  gives 
acceptable  results  has  yet  to  be  developed. 
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